Vec3f viewport_upper_left =
        center - m_focal_distance * w - viewport_u * 0.5f - viewport_v * 0.5f;

    pixel00_loc =
        viewport_upper_left + 0.5f * pixel_delta_u + 0.5f * pixel_delta_v;
    // Feature 1: Defocus Blur

    // Calculate the camera defocus disk basis vectors.
    auto defocus_radius = focus_dist * std::tan(Spherical::deg2rad(defocus_angle / 2));

   // FEATURE 2: APERATURE SHAPE SAMPLING 
inline Vec2f sample_ngon_aperture(int n) {
    // Decompose n-gon into n triangles from center
    // Task 1: Pick a random triangle (uniform)
    // Task 2: Sample uniformly within that triangle
    
    // Pick random triangle index
    int triangle_idx = std::min(int(randf() * n), n - 1);
    
    // Calculate vertices of the chosen triangle
    Vec2f v0(0.0f, 0.0f);  // Center
    
    // Two vertices on the perimeter
    float angle1 = 2.0f * M_PI * triangle_idx / n;
    float angle2 = 2.0f * M_PI * (triangle_idx + 1) / n;
    Vec2f v1(std::cos(angle1), std::sin(angle1));
    Vec2f v2(std::cos(angle2), std::sin(angle2));
    
    // Sample uniformly within this triangle using barycentric coordinates
    Vec2f b = sample_triangle(Vec2f(randf(), randf()));
    float alpha = b.x;
    float beta = b.y;
    float gamma = 1.0f - alpha - beta;
    
    // Return the sampled point
    return v0 * alpha + v1 * beta + v2 * gamma;
}

    /// Construct a new ray
    Ray() : o(), d(), mint(epsilon), maxt(infinity), time(T(0)), medium(nullptr) {}

    /// Construct a new ray
    Ray(const Vec &o, 
    const Vec &d, 
    T mint = Ray::epsilon, 
    T maxt = Ray::infinity, 
    T time = T(0),
    Medium* medium = nullptr) :
    o(o), d(d), mint(mint), maxt(maxt), time(time), medium(medium) {}

   parallel_for(
        blocked_range(0,              // begin
                                image.height(), // end
                                1               // block_size (process 1 row at a time)
                                ),

        // Lambda function that will be called for each block of rows
        [&](uint32_t begin, uint32_t end, int unit_index, int thread_index)
        {
            for (uint32_t y = begin; y != end; ++y)
            {
                ...

   bool sample_free_flight(const Ray3f &ray, int channel, Sampler &sampler, HitRecord &hit, Color3f &f,
                            Color3f &p) const override
    {
        // Get medium coefficients
        Color3f sigma_a, sigma_s, sigma_n;
        std::tie(sigma_a, sigma_s, sigma_n) = coeffs(ray.o);
        Color3f sigma_t = sigma_a + sigma_s + sigma_n;

        float t = -std::log(1 - sampler.next1f()) / sigma_t[channel];

        bool sampledMedium = t < ray.maxt;
        t = std::min(t, ray.maxt);

        // Compute transmittance to the sampled point
        Color3f Tr = exp(-sigma_t * t);

        // Compute PDF for the sampled distance
        Color3f density = sampledMedium ? (sigma_t * Tr) : Tr;
        float mixture_pdf = (density[0] + density[1] + density[2]) / 3.f;

        hit.t = t;
        hit.p = ray(t);

        if (sampledMedium) {
            // Scattering event in medium
            f *= Tr * sigma_s;
            p *= mixture_pdf;
            return true;
        } else {
            // Ray exited medium without scattering
            f *= Tr;
            p *= mixture_pdf;
            return false;
        }
    }

    Color3f total_transmittance(const Ray3f &ray, Sampler &sampler) const override
    {
        Color3f sigma_a, sigma_s, sigma_n;
        std::tie(sigma_a, sigma_s, sigma_n) = coeffs(ray.o);
        Color3f sigma_t = sigma_a + sigma_s + sigma_n;

        return exp(-sigma_t * ray.maxt);
    }
};

    vbool NanoVDBMedium::sample_free_flight(const Ray3f &ray, int channel, Sampler &sampler, HitRecord &hit, Color3f &f,
                                       Color3f &p) const
{
    // transform ray into volume's space
    Ray3f r = m_xform.inverse().ray(ray);

    // convert to a NanoVDB ray, and check intersection
    auto world_ray =
        nanovdb::Ray(nanovdb::Vec3f(r.o.x, r.o.y, r.o.z), nanovdb::Vec3f(r.d.x, r.d.y, r.d.z), r.mint, r.maxt);

    // clip grid-space ray to bounding box
    if (!world_ray.clip(m_density_grid->worldBBox()))
    {
        hit.t = ray.maxt;
        hit.p = ray(hit.t);
        return false; // ray does not hit any of the medium
    }

    else
    {
        // Delta tracking / null-scattering implementation for heterogeneous medium
        float t = world_ray.t0();

        while (true) {
            // Sample free-flight distance using majorant (m_total)
            t -= std::log(1 - sampler.next1f()) / m_total[channel];

            // Check if we've exited the volume
            if (t >= world_ray.t1())
            {
                // Ray exited medium without scattering
                hit.t = ray.maxt;
                hit.p = ray(hit.t);
                // Transmittance is accumulated through null-scattering in the integrator
                // Don't modify f and p - they track the cumulative weight and pdf
                return false;
            }

            // Evaluate coefficients at sampled position
            Vec3f pos = ray(t);
            Color3f sigma_a, sigma_s, sigma_n;
            std::tie(sigma_a, sigma_s, sigma_n) = coeffs(pos);

            // Compute probabilities for real vs null interaction
            float sigma_real = sigma_a[channel] + sigma_s[channel];
            float p_real = sigma_real / m_total[channel];

            // Russian roulette to determine interaction type
            if (sampler.next1f() < p_real)
            {
                float p_scatter = sigma_s[channel] / sigma_real;

                if (sampler.next1f() < p_scatter)
                {
                    hit.t = t;
                    hit.p = pos;
                    hit.mat = nullptr;

                    // Path weight and PDF for delta tracking
                    // The sigma_s cancels with the pdf in the MIS computation
                    f *= sigma_s;
                    p *= m_total;

                    return true;
                }
                else
                {
                    // Absorption event - path terminates
                    f *= Color3f(0.f);
                    return false;
                }
            }
            // else: null-scattering event - continue sampling
        }
    }
}

Color3f NanoVDBMedium::total_transmittance(const Ray3f &ray, Sampler &sampler) const
{
    // DONE: implement this function if you want to support nanovdb media. This can be identical to
    // HomogeneousMedium::total_transmittance()
    Ray3f r = m_xform.inverse().ray(ray);
    auto world_ray = nanovdb::Ray(nanovdb::Vec3f(r.o.x, r.o.y, r.o.z), nanovdb::Vec3f(r.d.x, r.d.y, r.d.z), r.mint, r.maxt);
    if (!world_ray.clip(m_density_grid->worldBBox()))
    {
        return Color3f(1.0f); // ray does not hit any of the medium
    }

    Color3f transmittance(1.0f);
    float t = world_ray.t0();

    while (true) {
        float mu_t = std::max({m_total.x, m_total.y, m_total.z});
        t -= std::log(1.0f - sampler.next1f()) / mu_t;

        if (t >= world_ray.t1())
            break;

        Vec3f pos = ray(t);
        Color3f sigma_a, sigma_s, sigma_n;
        std::tie(sigma_a, sigma_s, sigma_n) = coeffs(pos);

        transmittance *= sigma_n / m_total;
    }

    return transmittance;
}